Folding Kinetics of Proteins and Cold Denaturation. 
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Folding kinetics of a lattice model of protein is studied. It uses the Random Energy Model for the 
intrachain couplings and a temperature dependent free energy of solvation derived from a realistic 
hydration model of apolar solutes. The folding times are computed using Monte Carlo simulations 
in the region of the phase diagram where the chain occurs in the native structure. These folding 
times are roughly equals for the temperatures of cold and warm denaturation for a large range 
of solvent quality. Between these temperatures, the folding times reach maxima and thus, at low 
temperatures, the kinetics of the chain always speeds up as the temperature is decreased. The study 
of the conformational space as function of the temperature permits to elucidate this phenomenon. 
At low temperature, it shows that the activation barriers of the system decrease faster than the 
temperature as the temperature is decreased. At high temperature, the rate of the barriers over the 
temperature decreases as the temperature is increased because the height of the barrier is almost 
constant. 

PACS numbers: 87.14.et 



Proteins are very long molecular chains built with 
given sequences of amino- acids. Under biological solvent 
conditions, a protein occurs in a unique, native, compact 
form [Ij] and an important goal of theoretical physics is 
to understand how a chain finds its native structure in a 
reasonable biological time. 

Lattice models, in which the amino acids of the chain 
are located on the vertices of a two or three-dimensional 
lattice, are widely used to study protein folding. In 
the Random Energy Model (REM) [3, i, i, i, ij, the 
couplings between monomers are chosen at random in 
a Gaussian distribution centered on a negative value. 
It leads to an energy spectrum where a few (all com- 
pact) conformations lye in the bottom discrete part of 
the spectrum while the large majority of the conforma- 
tions belongs to a quasi-continuous top part. This spec- 
trum may well mimics that of proteins. However, al- 
though REM explains some features of the proteins, it 
is independent on the temperature and it fails to re- 
produce a general feature specific to proteins: the cold 
denaturation0, i, i, [H . 

The cold denaturation has been first observed by Pri- 
valov [TTI for Myoglobin which is in its native form be- 
tween a temperature of warm denaturation T^, and a 
temperature of cold denaturation T^. Above or be- 
low Tc, the protein is in a denatured state where a lot of 
conformations are relevant. These temperatures are very 
sensitive to the pH of the solvent. 

Under physiological conditions, the proteins strongly 
interact with the solvent [H, IH, H, [H, [H, [13 and then 
any simulation of protein folding must consider a realis- 
tic solvent effect on the chain conformations. Recently 
[isl [Tgl , the temperature dependence of the hydrophobic 
effect has been introduced in the couplings of REM us- 
ing realistic solvent model based on a qualitative study 
of the energy spectra of the pure solvent and of the sol- 
vent around a monomer [20j . As a result, the first phase 



diagram of protein, where both warm and cold denat- 
urations occur has been calculated[l^ showing a very 
good accordance with experimental data. On the other 
hand, several works has been published on the subject 
providing alternative models which were able to exhibits 
theoretically the cold denaturation [2l|, [H [H, [13, [H [111 
In this paper, this protein model is first reminded. The 
different contributions to the effective couplings between 
monomers are shown as functions of the temperature. 
This result gives an insight of the modification of confor- 
mational space as the temperature is varied. Thus, the 
kinetic properties of this model are studied in the native 
region. Folding times are computed versus a solvation 
parameter and the temperature, using Monte Carlo sim- 
ulations 127|, i2^]. Moreover, the native conformation, the 
kinetic trap and the transition state are determined by 
a study of the phase space. Last, the unusual behavior 
of the folding times at low temperature is elucidated by 
taking into account the temperature dependence of the 
free energy of the activation barrier. 



I. PROTEIN MODEL 

The chain is represented by a string of A'' beads, (here 
N — 16), located on a square two-dimensional lattice. 
For a given chain conformation, each empty lattice ver- 
tice is considered to be a solvent site. It is filled up 
with four solvent cells pointing towards the four direc- 
tions (see fig HI). Thus, a solvent cell interacts either with 
a monomer or with another solvent cell and its thermo- 
dynamic properties are determined by the type of this 
interaction. In such a model, the number of solvent sites 
and then the number of solvent cells are constant be- 
cause the chain length is fixed. Hence, the volume of 
the solvent does not depend on the chain conformation. 
Moreover, here, a unique parameter gives an insight 
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of the solvent quality. 
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II. EFFECTIVE HAMILTONIAN 

The system is composed by the chain and a very large 
number of solvent cells. It is at equilibrium with a bath 
at temperature T and the probability of occurrence of a 
chain conformation m in interaction with the solvent is : 



The effective Hamiltonian of conformation m takes into 
account of all the lattice links and has the following form : 

JV N 

i>j + l i=l 



water-water 



intrachain 



solvation 



where Bij is the coupling parameter between the 
monomers i and j, chosen at random in a Gaussian dis- 
tribution centered on with a standard deviation equal 
to 2 and A^™'' equals 1 if the monomers i and j are first 
neighbors on the lattice and otherwise. fi{T) is the free 
energy of a solvent cell in interaction with monomer i and 
fs{Bs,T) is the free energy of a pure solvent cell, tt.^™'' is 
the number of links between monomer i and solvent nods 
and ni™'' is the number of solvent vertices bonds. 

The expression of the constant number of total lattice 
links (equals to Ki) and the fact that each monomer 
always creates 4 links leads to write the two conservation 
equations of the model : 
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One deduces easily : 
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Owing to these relations, and after removing the constant 
terms, the effective Hamiltonian is rewritten as : 
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It takes the form of the usual Hamiltonian without sol- 
vent effect with effective coupling parameters which now 
depend on the solvent properties {Bg) and the tempera- 
ture : 

Btf{B,, T) = - /.(T) - /, (T) + 2/.(B„ T) 

The model used in this work is that introduced in 
ref.l9. The main difference with that used in ref.18 comes 
from the factor 2 associated to the free energy of a pure 
solvent cell in the form of the effective couplings. This 
function guarantees that the total number of solvent cells 
(i.e. the volume of the solvent) is a constant whatever 
the chain conformation. In other words, the creation of 
an link between monomers and j involves the removing of 
an interaction between residues i in one hand and j in an 
another hand with a solvent cell each. These two solvent 
cells becomes pure solvent. 

Obviously, the final form of the effective couplings de- 
pends on the solvent model used to calculated fs and 



III. SOLVATION MODEL 

The solvation free energy calculation are based on re- 
sults of a study of the hydrophobic effect undertaken 
by Dill and coworkersj20j . They used the Mercedes 
Benz model of water |29j and a simple adaptation of 
the two-states model of Mullerjs^l extended by Lee and 



Graziano 3l| to give a physical picture of the hydropho- 
bic effect in terms of two energy spectra. The first one 
is associated to the pure solvent and the other one to 
water molecules in interaction with an apolar solute. In 
the spectrum associated to pure solvent, the energy gap 
between the ground and the low lying exited states is 
small. Here, these two close energy states are gathered 
in a unique one of energy Bg, A^"-folds degenerate (see 
fig. It comes : 



fs{Bs,T) = Bs-aT\nNs 



with 



a <1 



In the other spectrum, the energy gap is larger and 
the exited state is more degenerate. This picture is re- 
duced further in the spirit of the REM and the spec- 
trum is spread out (see figlUb)). For each monomer i, 
A'^s values of the solvation energies are draw at random 
in a Gaussian distribution centered on with standard 
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FIG. 2: .(a) Energy spectrum of a cell of pure solvent with a 
unique A^^ fold degenerated level of energy Bs- (b) Energy 
spectrum of a cell of solvent in interaction with a monomer. 
Ns energy values are draw at random in a Gaussian distribu- 
tion centered on 0. 5™'° is specific to each monomer. 



deviation equals 2 and the minimum value of each set 
of solvation energies, specific to each monomer is deter- 
mined and noted i?™'". The free energy of solvent cell in 
interaction with the monomer i is then : 



-Tin 



giB)eM~B/T)dB 



where g{B) is the density of energy states, i.e., a Gaussian 
truncated at B™™. Thus, one has : 




dB 



IV. HYDRATION RESULTS 

The solvation parameters used in this work are a = 
a = 2 and Ns = 10^. 

As it has already pointed out in a another work[l£ 
choosing such solvent parameters, one has : 
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Then, using a — 0.5 in the free energy of the pure solvent 
is equivalent to the previous work[l8j where the solvent 
parameter has only been divided by 2. 

To well understand the kinetic of folding of the chain 
present later, we first focus on the different contributions 
to the effective couplings. Figure [3] shows that for large 
values of Bg , the effective couplings are repulsive, what- 
ever the temperature, because Bij + 2fs > fi + fj- Thus, 
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FIG. 3: Curves of the different contributions to the effective 
coupling between the monomers i=l and j=4 as function of 
the temperature for several values of the solvent quality. The 
two temperatures for which the coupling vanishes are shown 
for Bs = -6.0. 



this condition corresponds to a good solvent in the sense 
that the monomers are preferably exposed to the solvent. 
For small values oiBg, the couplings are attractive at low 
temperature, modelling bad solvent condition, and then 
repulsive at high temperature. Last, for intermediate val- 
ues of Bg, the curves of Bij + 2fs intersect those of fi + fj 
twice at two temperatures noted T-7 and T^. The effec- 
five couplings are repulsive for T < or T > and 

attractive for T,7 < T < . 

''J ''J 

These results show that the free energies of transfer 
of the residues into water (here Sfi = fi — fs) present 
maxima for temperatures, depending on Bs, between 
and This is in good agreement with studies of the 
temperature dependence of the hydrophobic interaction 
in protein folding. As an example, first, Baldwinjs^l 
showed from calorimetric data, that the transfer of hy- 
drocarbons in water always exhibits a temperature (de- 
noted as Ts) for which the entropy of transfer reaches 
zero {AS{Ts) = 0). Using the fundamental thermody- 
namical relation AS* = —dAF/dT, it is clear that the 
free energy of transfer, A_F of the hydrocarbons reaches 
a maximum at T^. 



PROTEIN THERMODYNAMICS 



It must be noted that, for given solvent quality and 
temperature, each intrachain couplings, Bij, are differ- 
ent from each other. The functions fi{T) are also spe- 
cific to each residue because the i?™'" depend on the 
monomers. Thus, some couplings are more attractive 
than other ones. As a results, the ground state of the ef- 
fective Hamiltonian spectra are non degenerated when 
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the chain is in interaction with a bad quahty solvent 
{Bs ^ 0) . It corresponds to the native, more maxi- 
maUy compact, conformation of the sequence. The native 
structure (Nat) is determined by a full enumeration of the 
conformational space of the chain in interaction with a 
bad solvent. For several values of Bs and T, the proba- 
bility of occurrence of Nat is calculated. If P^^*" > 1/2, 
the chain is in the native phase. 




FIG. 4: Logarithm of the folding time of the chain versus 
Bs and T in the native region for —7.0 < Bs < —5.75 and 
n{Bs) <T < T^{Bs). 



VI. PROTEIN KINETICS 

The folding time, ifoid, only defined in the native phase, 
is the mean Monte Carlo steps needed to reach Nat for 
the first time, averaged over 1000 trajectoriesjs^ . Each 
trajectory starts with a random conformation and Monte 
Carlo simulations using the corner flip, the tail and the 
crankshaft moves used in [il, Is^l, are performed. The 
folding times are plotted versus B^ and T in figlH 

For -5.75 > J5,, > -6.4, one has ifoid(Tc) « ifoid(T^)- 
That is to say, the folding times are always the same at 
the temperatures of denaturation in this case (see figE]). 
This property may be understood as follow. By defini- 
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The inset of figH] shows the phase diagram of this se- 
quence. For Bg > —5.75, the chain is always in the de- 
natured phase due to the repulsive couplings and then a 
lot of different structures are relevant. For —7.0 < Bg < 
—5.75, the chain occurs in the native state between Tc 
and Tw depending on Bs- Below Tc and above Tw, the 
chain is denatured. Obviously the cold denaturation, due 
to the change in the sign of the effective couplings of Nat, 
occurs for Tc ~ T^ . Moreover, with the parameters used 

in this work, one also has T^ « T^ . For temperature 
smaller than Tc, the couplings are mainly repulsive and 
then only the numerous chain structures without con- 
tact are relevant. Thus, in the cold phase, the probabil- 
ity of occurrence of the Native structure becomes very 
small. This phase is well denatured because, if it is a 
glassy phase, the equilibrium probability of occurrence 
of the native structure would be rather large (and kinet- 
ics would be very slow). In the cold denatured phase, 
the set of the extended structures becomes the relevent 
sampling of conformations. Hence, kinetics would not 
converge towards the Native structure but would diffuse 
freely among the extended structures subset. 

At Bg = —7.0, the couplings are always attractive at 
low temperature. The cold denaturation disappears and 
for Bs < —7.0, only the warm transition remains. A 
critical point occurs for Bs'^'^ = —5.75 and T^'^^ = 0.53. 
At this point, one has T^ « T^. 



FIG. 5: . Folding times at the temperatures of denaturation 
as a function of the solvent quality. The filled circle are for 
the temperature of cold denaturation and the -I- are for the 
warm denaturation. 

tion, the equilibrium probability of occurrence of the na- 
tive structure equals 1/2 at Tc and T^j. Moreover, it may 
be seen on the example of figEl that, for Bs = —6.25, 
the temperatures of transition are T^ — 0.83 and Tc = 
0.24 and the effective Hamiltonian of the native struc- 
ture for these temperatures are TC^g^lTc) — —4.2 and 
n^i^T^) = -14.2 leading to 

HT{Tc)/Tc^nT{T^)/T^ 

As, 

eM~'H^i'{Tc)/Tc] _ exp[-H,y(r^)/r^] 
Z{Tc) Z{T^) 

one deduces the equality of the partition functions, 
Z{Tc) = Z{Tiu). Thus, the probability of occurrence of 
each extended structure, having a zero effective Hamil- 
tonian, is simply one over the partition function and 
is the same at Tc and T^,. Rationally, to satisfy to 
the equality of the partition functions, one supposes 
that the conformational effective Hamiltonians satisfy to 
ni'^\Tc)/Tc w nig\T^,)/T^ whatever the structure to. 
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The configurational spaces where the effective Hamilto- 
nians are roughly linearly scaled by the temperatures 
may thus be expected for these temperatures. This re- 
sult holds whatever the solvent quality. To be more pre- 
cise, the equilibrium probabilities of each structure re- 
mains almost constant at the temperatures of denatura- 
tion when the solvent quality is varied. As a consequence, 
the transition rates between two conformations m and n, 
wim n,T) ^ exp[(H^^)(T) - n'';^\T)) /2T], are al- 
most constant at the denaturations temperatures for all 
Bs > —6.4, leading to this (quasi) equality of the folding 
times at Tc and 

For Bs —7.0, the temperature of cold denaturation 
tends toward 0. The conformational landscape still ex- 
hibits a well pronounced effective Hamiltonian minimum 
for the native structure, but as the temperature becomes 
very small the transition rates becomes either quasi null 
or huge. It takes huge amount of time to overcome some 
local energetic barriers. Thus the folding time tends to 
infinity at Tc because the kinetics is then frozen. 

Then, folding time at Tc increases rapidly as Bs de- 
creases from -6.4 to -7.0. 

For each value of Bs, the folding time goes through a 
maximum value at a temperature noted T*{Bs). The 
region of the curves where T > T*, simply confirms 
that the time needed for a random conformation to reach 
the native structure increases as the temperature is de- 
creased. However, the behavior of the curves for T <T* 
is less usual. It corresponds to an increasing of the speed 
of the kinetics with respect to a decrease of the temper- 
ature. This result is not in accordance with standard 
Monte Carlo simulations and its explanation is given be- 
low. 

For each temperature, the kinetic trap [33] (shown in 
figll]) is found using Monte Carlo algorithm by noting 
the most occurrent conformation during the 1000 trajec- 
tories performed to calculated the folding time. It is ob- 
served that the trap do not depend on the temperature. 
The trap present a great similarity with the native struc- 
ture. Six contacts among the nine contacts of Nat are 
observed in trap. Thus, its effective Hamiltonian is very 
close from that of Nat showing that this structure could 
be a good candidate for the trap of the system. Second, 
there is no physical pathway from trap to Nat keeping 
unbroken the whole set of contacts. That is to say, in a 
typical trajectory from the trap to Nat, each contact of 
the trap has to be broken and some of them created again 
in an other global arrangement. Obviously, they are not 
simultaneously broken. If it were the case, the trajec- 
tories would contain some extended structures without 
contacts which would have a very great difference of ef- 
fective Hamiltonian with the trap. Consequently, even 
step by step, the transition rate for such sub-pathway 
would be quasi-null. Thus, in more favorable pathways, 
some contacts are broken while some others are created 
and the chain walks through local rearrangements toward 
Nat without reaching any extended conformation. A par- 
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FIG. 6: Inset: schematic free energy landscape of a folding 
pathway from the trap (Tp) to the native (Nat) structure 
through the transition state (TS) with the activation barrier 
(arrows) l^F^{T) = H^g^' -^^^^ Main: free energy of the 
native (square) , the trap (filled circle) conformations and the 
transition state (-1-) versus temperature for Bs = —6.25. The 
arrows show some free energy barriers A_F*(r). 

ticular structure of the sequence is the transition state. 
It is determined by performing a new simulation of 1000 
trajectories where the trap is always the first conforma- 
tion. For each trajectory, the conformation of highest 
value of the effective Hamiltonian is considered as a pos- 
sible transition state. The transition state (TS) is the 
structure with the lowest value of effective Hamiltonian 
among the sampling of possible transition states collected 
over the 1000 trajectories (see inset of figlB])- Among all 
the possible pathways from trap to Nat, those passing 
through TS are the less energetically costly. One must 
note that this is not the usual definition of the transition 
state adopted in the theory of protein folding where the 
TS is not a unique structure but an ensemble of config- 
urations of highest free energy along the path or paths 
between unfolded macrostates and the native structure 
[3^. Here, the defintion of the TS is that used in the 
theory of the simple gaz chemical reactions. Moreover, 
as explained above, it can be seen that TS has a weak 
similarity with Nat or with the trap. 

The key role of the trap in the folding process is ex- 
plained as follow. Some trajectories, starting by a ran- 
dom conformations, fall in the trap valley. To leave this 
structure the chain has to overcome the largest barrier of 
the system. On the other hand, the trajectories which do 
not reach the trap valley, walk down the native structure, 
passing quickly over smaller barriers. Thus, only the trap 
activation barrier is considered in the folding analysis. 

For Bs = —6.25, the values of the effective Hamil- 
tonian for Nat, the trap and TS are reported on figEl 
The activation barrier between the trap and Nat struc- 
tures, AF^ = H^g^'' — Ti''^^'^^\ shown on figUl obviously 
depends on the temperature. For T > 0.60, AF^ is al- 
most constant whereas it decreases very quickly with the 
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temperature for T < 0.45. Moreover, the kinetic theory 
indicates that the time needed to escape from the trap 
vaUey is in proportion to exp(AF-'-/T). Figure [7] showed 
the very sharp decreasing of AF^ /T with respect to a de- 
crease of the temperature for T > 0.45. One must noted 
the quasi equahty of this last quantity at the tempera- 
tures of denaturation. 



is decreased because AF^ decreases faster than T. For 
T > Ti, the rate AF^/T decreases as the temperature is 
increased because AF* is almost constant and, obviously, 
1/T decreases. As, the folding time is in proportion to 
exp(AF-'-(r)/T), this result permits to elucidate why the 
lower the temperature, the faster is the kinetics at low 
temperature and the higher the temperature, the faster 
also is the kinetics at high temperature. 



--AF^/T 




FIG. 7: Logarithm of the folding time (filled square) and 
AF^T)/T (filled circle) versus the temperature for Bs = 
-6.25. 

The values of AF^ /T and that of the folding time 
starting from a random structure, reported in figEl 
present maxima at the temperature, Ti « 0.45. For 
T < Ti, the rate AF^/T decreases as the temperature 



VII. CONCLUSION 

Our model predicts and explains unusual behavior of 
the folding times due to the flatness of the conformation 
space as the temperature decreases towards Tc- In the 
native region, far from the denaturation temperatures, 
P^^^ tends to 1, leading to a very rough conformational 
landscape. On the opposite, around Tc, effective Hamil- 
tonian of all conformations are quasi equal and the free 
energy landscape is quasi flat. For low temperatures, the 
difference of effective Hamiltonian between two confor- 
mations decreases faster than the temperature. This lev- 
eling of the conformational space close to the cold transi- 
tion leads to a drastic decrease of the free energy barriers 
of the configurational space leading to an acceleration of 
the kinetics as the temperature decreases. Last and un- 
fortunately, at the present time, they are no experimen- 
tal results available for the folding kinetics of proteins at 
temperatures close to that of cold denaturation. which 
should be compared to this theoretical prediction. 
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